## Description #############################################################################
#
# Tests related to ECI to ECI transformations.
#
## References ##############################################################################
#
# [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications. Microcosm Press,
#     Hawthorn, CA, USA.
#
############################################################################################

# Get the current EOP.
#
# TODO: The EOP obtained from IERS website does not match the values in the examples in [1].
# However, this should be enough, because 1) the individual functions at the low level are
# tested using the same values of [1], and 2) the difference is smaller than 30 cm.

eop_iau1980  = read_iers_eop("../eop_IAU1980.txt",  Val(:IAU1980))
eop_iau2000a = read_iers_eop("../eop_IAU2000A.txt", Val(:IAU2000A))

# == File: ./src/transformations/eci_to_eci.jl =============================================

############################################################################################
#                                       IAU-76 / FK5                                       #
############################################################################################

# -- Functions: r_eci_to_eci ---------------------------------------------------------------

# == GCRF <=> J2000 ========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895790   i + 6123.01140070   j + 6378.13692820   k [km]
#
# one gets:
#
#   r_j2000 = 5102.50960000   i + 6123.01152000   j + 6378.13630000   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci GCRF <=> J2000" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == GCRF => J2000 =====================================================================

    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]

    # -- DCM -------------------------------------------------------------------------------

    D_j2000_gcrf = r_eci_to_eci(GCRF(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = D_j2000_gcrf * r_gcrf

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_j2000_gcrf = r_eci_to_eci(Quaternion, GCRF(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = vect(q_j2000_gcrf \ r_gcrf * q_j2000_gcrf)

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-4

    # == J2000 => GCRF =====================================================================

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]

    # -- DCM -------------------------------------------------------------------------------

    D_GCRF_J2000 = r_eci_to_eci(J2000(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = D_GCRF_J2000 * r_j2000

    @test r_gcrf[1] ≈ +5102.50895790 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_gcrf_j2000 = r_eci_to_eci(Quaternion, J2000(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = vect(q_gcrf_j2000 \ r_j2000 * q_gcrf_j2000)

    @test r_gcrf[1] ≈ +5102.50895790 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol = 1e-4
end

# == GCRF <=> MOD ==========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895790   i + 6123.01140070   j + 6378.13692820   k [km]
#
# one gets:
#
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci GCRF <=> MOD" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == GCRF => MOD =======================================================================

    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]

    # -- DCM -------------------------------------------------------------------------------

    D_mod_gcrf = r_eci_to_eci(GCRF(), MOD(), JD_UTC, eop_iau1980)

    r_mod = D_mod_gcrf * r_gcrf

    @test r_mod[1] ≈ +5094.02837450 atol = 1e-4
    @test r_mod[2] ≈ +6127.87081640 atol = 1e-4
    @test r_mod[3] ≈ +6380.24851640 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_mod_gcrf = r_eci_to_eci(Quaternion, GCRF(), MOD(), JD_UTC, eop_iau1980)

    r_mod = vect(q_mod_gcrf \ r_gcrf * q_mod_gcrf)

    @test r_mod[1] ≈ +5094.02837450 atol = 1e-4
    @test r_mod[2] ≈ +6127.87081640 atol = 1e-4
    @test r_mod[3] ≈ +6380.24851640 atol = 1e-4

    # == MOD => GCRF =======================================================================

    r_mod = [5094.02837450; 6127.87081640; 6380.24851640]
    v_mod = [-4.7462630520; 0.7860140450; 5.5317905620]

    # -- DCM -------------------------------------------------------------------------------

    D_gcrf_mod = r_eci_to_eci(MOD(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = D_gcrf_mod * r_mod

    @test r_gcrf[1] ≈ +5102.50895790 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_gcrf_mod = r_eci_to_eci(Quaternion, MOD(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = vect(q_gcrf_mod \ r_mod * q_gcrf_mod)

    @test r_gcrf[1] ≈ +5102.50895790 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol = 1e-4
end

# == GCRF <=> TOD ==========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895790   i + 6123.01140070   j + 6378.13692820   k [km]
#
# one gets:
#
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci GCRF <=> TOD" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == GCRF => TOD =======================================================================

    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]

    # -- DCM -------------------------------------------------------------------------------

    D_tod_gcrf = r_eci_to_eci(GCRF(), TOD(), JD_UTC, eop_iau1980)

    r_tod = D_tod_gcrf * r_gcrf

    @test r_tod[1] ≈ +5094.51620300 atol = 1e-4
    @test r_tod[2] ≈ +6127.36527840 atol = 1e-4
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_tod_gcrf = r_eci_to_eci(Quaternion, GCRF(), TOD(), JD_UTC, eop_iau1980)

    r_tod = vect(q_tod_gcrf \ r_gcrf * q_tod_gcrf)

    @test r_tod[1] ≈ +5094.51620300 atol = 1e-4
    @test r_tod[2] ≈ +6127.36527840 atol = 1e-4
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-4

    # == TOD => GCRF =======================================================================

    r_tod = [5094.51620300; 6127.36527840; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_gcrf_tod = r_eci_to_eci(TOD(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = D_gcrf_tod * r_tod

    @test r_gcrf[1] ≈ +5102.50895790 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_gcrf_tod = r_eci_to_eci(Quaternion, TOD(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = vect(q_gcrf_tod \ r_tod * q_gcrf_tod)

    @test r_gcrf[1] ≈ +5102.50895790 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol = 1e-4
end

# == GCRF <=> TEME =========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895790   i + 6123.01140070   j + 6378.13692820   k [km]
#
# one gets:
#
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci GCRF <=> TEME" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == GCRF => TEME ======================================================================

    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]

    # -- DCM -------------------------------------------------------------------------------

    D_teme_gcrf = r_eci_to_eci(GCRF(), TEME(), JD_UTC, eop_iau1980)

    r_teme = D_teme_gcrf * r_gcrf

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-4
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-4
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_teme_gcrf = r_eci_to_eci(Quaternion, GCRF(), TEME(), JD_UTC, eop_iau1980)

    r_teme = vect(q_teme_gcrf \ r_gcrf * q_teme_gcrf)

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-4
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-4
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-4

    # == TEME => GCRF ======================================================================

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_gcrf_teme = r_eci_to_eci(TEME(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = D_gcrf_teme * r_teme

    @test r_gcrf[1] ≈ +5102.50895790 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_gcrf_teme = r_eci_to_eci(Quaternion, TEME(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = vect(q_gcrf_teme \ r_teme * q_gcrf_teme)

    @test r_gcrf[1] ≈ +5102.50895790 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol = 1e-4
end

# == J2000 <=> MOD =========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC     = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#
# one gets:
#
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#
# If not EOP correction is used, then we get:
#
#   r_mod = 5094.02901670   i + 6127.87093630   j + 6380.24788850   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci J2000 <=> MOD" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == J2000 => MOD ======================================================================

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]

    # -- DCM -------------------------------------------------------------------------------

    D_MOD_J2000 = r_eci_to_eci(J2000(), MOD(), JD_UTC, eop_iau1980)

    r_mod = D_MOD_J2000 * r_j2000

    @test r_mod[1] ≈ +5094.02837450 atol = 1e-4
    @test r_mod[2] ≈ +6127.87081640 atol = 1e-4
    @test r_mod[3] ≈ +6380.24851640 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_mod_j2000 = r_eci_to_eci(Quaternion, J2000(), MOD(), JD_UTC, eop_iau1980)

    r_mod = vect(q_mod_j2000 \ r_j2000 * q_mod_j2000)

    @test r_mod[1] ≈ +5094.02837450 atol = 1e-4
    @test r_mod[2] ≈ +6127.87081640 atol = 1e-4
    @test r_mod[3] ≈ +6380.24851640 atol = 1e-4

    # == MOD => J2000 ======================================================================

    r_mod = [5094.02837450; 6127.87081640; 6380.24851640]

    # -- DCM -------------------------------------------------------------------------------

    D_j2000_mod = r_eci_to_eci(MOD(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = D_j2000_mod * r_mod

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_j2000_mod = r_eci_to_eci(Quaternion, MOD(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = vect(q_j2000_mod \ r_mod * q_j2000_mod)

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-4

    # == No EOP Corrections ================================================================

    # == J2000 => MOD ======================================================================

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]

    # -- DCM -------------------------------------------------------------------------------

    D_MOD_J2000 = r_eci_to_eci(J2000(), MOD(), JD_UTC)

    r_mod = D_MOD_J2000 * r_j2000

    @test r_mod[1] ≈ +5094.02901670 atol = 1e-7
    @test r_mod[2] ≈ +6127.87093630 atol = 1e-7
    @test r_mod[3] ≈ +6380.24788850 atol = 1e-7

    # -- Quaternion ------------------------------------------------------------------------

    q_mod_j2000 = r_eci_to_eci(Quaternion, J2000(), MOD(), JD_UTC)

    r_mod = vect(q_mod_j2000 \ r_j2000 * q_mod_j2000)

    @test r_mod[1] ≈ +5094.02901670 atol = 1e-7
    @test r_mod[2] ≈ +6127.87093630 atol = 1e-7
    @test r_mod[3] ≈ +6380.24788850 atol = 1e-7

    # == MOD => J2000 ======================================================================

    r_mod = [5094.02901670; 6127.87093630; 6380.24788850]

    # -- DCM -------------------------------------------------------------------------------

    D_j2000_mod = r_eci_to_eci(MOD(), J2000(), JD_UTC)

    r_j2000 = D_j2000_mod * r_mod

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-7

    # -- Quaternion ------------------------------------------------------------------------

    q_j2000_mod = r_eci_to_eci(Quaternion, MOD(), J2000(), JD_UTC)

    r_j2000 = vect(q_j2000_mod \ r_mod * q_j2000_mod)

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-7
end

# == J2000 <=> TOD =========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#
# one gets:
#
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#
# If not EOP correction is used, then we get:
#
#   r_tod = 5094.51478040   i + 6127.36646120   j + 6380.34453270   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci J2000 <=> TOD" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == J2000 => TOD ======================================================================

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]

    # -- DCM -------------------------------------------------------------------------------

    D_TOD_J2000 = r_eci_to_eci(J2000(), TOD(), JD_UTC, eop_iau1980)

    r_tod = D_TOD_J2000 * r_j2000

    @test r_tod[1] ≈ +5094.51620300 atol = 1e-4
    @test r_tod[2] ≈ +6127.36527840 atol = 1e-4
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_tod_j2000 = r_eci_to_eci(Quaternion, J2000(), TOD(), JD_UTC, eop_iau1980)

    r_tod = vect(q_tod_j2000 \ r_j2000 * q_tod_j2000)

    @test r_tod[1] ≈ +5094.51620300 atol = 1e-4
    @test r_tod[2] ≈ +6127.36527840 atol = 1e-4
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-4

    # == TOD => J2000 ======================================================================

    r_tod = [5094.51620300; 6127.36527840; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_j2000_tod = r_eci_to_eci(TOD(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = D_j2000_tod * r_tod

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_j2000_tod = r_eci_to_eci(Quaternion, TOD(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = vect(q_j2000_tod \ r_tod * q_j2000_tod)

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-4

    # == No EOP corrections ================================================================

    # == J2000 => TOD ======================================================================

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]

    # -- DCM -------------------------------------------------------------------------------

    D_TOD_J2000 = r_eci_to_eci(J2000(), TOD(), JD_UTC)

    r_tod = D_TOD_J2000 * r_j2000

    @test r_tod[1] ≈ +5094.51478040 atol = 1e-7
    @test r_tod[2] ≈ +6127.36646120 atol = 1e-7
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-7

    # -- Quaternion ------------------------------------------------------------------------

    q_tod_j2000 = r_eci_to_eci(Quaternion, J2000(), TOD(), JD_UTC)

    r_tod = vect(q_tod_j2000 \ r_j2000 * q_tod_j2000)

    @test r_tod[1] ≈ +5094.51478040 atol = 1e-7
    @test r_tod[2] ≈ +6127.36646120 atol = 1e-7
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-7

    # == TOD => J2000 ======================================================================

    r_tod = [5094.51478040; 6127.36646120; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_j2000_tod = r_eci_to_eci(TOD(), J2000(), JD_UTC)

    r_j2000 = D_j2000_tod * r_tod

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-7

    # -- Quaternion ------------------------------------------------------------------------

    q_j2000_tod = r_eci_to_eci(Quaternion, TOD(), J2000(), JD_UTC)

    r_j2000 = vect(q_j2000_tod \ r_tod * q_j2000_tod)

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-7
end

# == J2000 <=> TEME ========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#
# one gets:
#
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci J2000 <=> TEME" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == J2000 => TEME =====================================================================

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]

    # -- DCM -------------------------------------------------------------------------------

    D_TEME_J2000 = r_eci_to_eci(J2000(), TEME(), JD_UTC, eop_iau1980)

    r_teme = D_TEME_J2000 * r_j2000

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-4
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-4
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_teme_j2000 = r_eci_to_eci(Quaternion, J2000(), TEME(), JD_UTC)

    r_teme = vect(q_teme_j2000 \ r_j2000 * q_teme_j2000)

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-4
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-4
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-4

    # == TEME => J2000 =====================================================================

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_j2000_teme = r_eci_to_eci(TEME(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = D_j2000_teme * r_teme

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_j2000_teme = r_eci_to_eci(Quaternion, TEME(), J2000(), JD_UTC)

    r_j2000 = vect(q_j2000_teme \ r_teme * q_j2000_teme)

    @test r_j2000[1] ≈ +5102.50960000 atol = 1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol = 1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol = 1e-4
end

# == MOD <=> TOD ===========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#
# one gets:
#
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#
# If not EOP correction is used, then we get:
#
#   r_mod = 5094.02901670   i + 6127.87093630   j + 6380.24788850   k [km]
#
#   r_tod = 5094.51478040   i + 6127.36646120   j + 6380.34453270   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci MOD <=> TOD" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == MOD => TOD ========================================================================

    r_mod = [5094.02837450; 6127.87081640; 6380.24851640]

    # -- DCM -------------------------------------------------------------------------------

    D_tod_mod = r_eci_to_eci(MOD(), JD_UTC, TOD(), JD_UTC, eop_iau1980)

    r_tod = D_tod_mod * r_mod

    @test r_tod[1] ≈ +5094.51620300 atol = 1e-4
    @test r_tod[2] ≈ +6127.36527840 atol = 1e-4
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_tod_mod = r_eci_to_eci(Quaternion, MOD(), JD_UTC, TOD(), JD_UTC, eop_iau1980)

    r_tod = vect(q_tod_mod \ r_mod * q_tod_mod)

    @test r_tod[1] ≈ +5094.51620300 atol = 1e-4
    @test r_tod[2] ≈ +6127.36527840 atol = 1e-4
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-4

    # == TOD => MOD ========================================================================

    r_tod = [5094.51620300; 6127.36527840; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_mod_tod = r_eci_to_eci(TOD(), JD_UTC, MOD(), JD_UTC, eop_iau1980)

    r_mod = D_mod_tod * r_tod

    @test r_mod[1] ≈ +5094.02837450 atol = 1e-4
    @test r_mod[2] ≈ +6127.87081640 atol = 1e-4
    @test r_mod[3] ≈ +6380.24851640 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_mod_tod = r_eci_to_eci(Quaternion, TOD(), JD_UTC, MOD(), JD_UTC, eop_iau1980)

    r_mod = vect(q_mod_tod \ r_tod * q_mod_tod)

    @test r_mod[1] ≈ +5094.02837450 atol = 1e-4
    @test r_mod[2] ≈ +6127.87081640 atol = 1e-4
    @test r_mod[3] ≈ +6380.24851640 atol = 1e-4

    # == No EOP corrections ================================================================

    # == MOD => TOD ========================================================================

    r_mod = [5094.02901670; 6127.87093630; 6380.24788850]

    # -- DCM -------------------------------------------------------------------------------

    D_tod_mod = r_eci_to_eci(MOD(), JD_UTC, TOD(), JD_UTC)

    r_tod = D_tod_mod * r_mod

    @test r_tod[1] ≈ +5094.51478040 atol = 1e-7
    @test r_tod[2] ≈ +6127.36646120 atol = 1e-7
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-7

    # -- Quaternion ------------------------------------------------------------------------

    q_tod_mod = r_eci_to_eci(Quaternion, MOD(), JD_UTC, TOD(), JD_UTC)

    r_tod = vect(q_tod_mod \ r_mod * q_tod_mod)

    @test r_tod[1] ≈ +5094.51478040 atol = 1e-7
    @test r_tod[2] ≈ +6127.36646120 atol = 1e-7
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-7

    # == TOD => MOD ========================================================================

    r_tod = [5094.51478040; 6127.36646120; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_mod_tod = r_eci_to_eci(TOD(), JD_UTC, MOD(), JD_UTC)

    r_mod = D_mod_tod * r_tod

    @test r_mod[1] ≈ +5094.02901670 atol = 1e-7
    @test r_mod[2] ≈ +6127.87093630 atol = 1e-7
    @test r_mod[3] ≈ +6380.24788850 atol = 1e-7

    # -- Quaternion ------------------------------------------------------------------------

    q_mod_tod = r_eci_to_eci(Quaternion, TOD(), JD_UTC, MOD(), JD_UTC)

    r_mod = vect(q_mod_tod \ r_tod * q_mod_tod)

    @test r_mod[1] ≈ +5094.02901670 atol = 1e-7
    @test r_mod[2] ≈ +6127.87093630 atol = 1e-7
    @test r_mod[3] ≈ +6380.24788850 atol = 1e-7
end

# == MOD <=> TEME ==========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#
# one gets:
#
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#
# If not EOP correction is used, then we get the same result by using:
#
#   r_mod = 5094.02901670   i + 6127.87093630   j + 6380.24788850   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci MOD <=> TEME" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == MOD => TEME =======================================================================

    r_mod = [5094.02837450; 6127.87081640; 6380.24851640]

    # -- DCM -------------------------------------------------------------------------------

    D_teme_mod = r_eci_to_eci(MOD(), JD_UTC, TEME(), JD_UTC, eop_iau1980)

    r_teme = D_teme_mod * r_mod

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-4
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-4
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_teme_mod = r_eci_to_eci(Quaternion, MOD(), JD_UTC, TEME(), JD_UTC, eop_iau1980)

    r_teme = vect(q_teme_mod \ r_mod * q_teme_mod)

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-4
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-4
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-4

    # == TEME => MOD =======================================================================

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_mod_teme = r_eci_to_eci(TEME(), JD_UTC, MOD(), JD_UTC, eop_iau1980)

    r_mod = D_mod_teme * r_teme

    @test r_mod[1] ≈ +5094.02837450 atol = 1e-4
    @test r_mod[2] ≈ +6127.87081640 atol = 1e-4
    @test r_mod[3] ≈ +6380.24851640 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_mod_teme = r_eci_to_eci(Quaternion, TEME(), JD_UTC, MOD(), JD_UTC, eop_iau1980)

    r_mod = vect(q_mod_teme \ r_teme * q_mod_teme)

    @test r_mod[1] ≈ +5094.02837450 atol = 1e-4
    @test r_mod[2] ≈ +6127.87081640 atol = 1e-4
    @test r_mod[3] ≈ +6380.24851640 atol = 1e-4

    # == No EOP corrections ================================================================

    # == MOD => TEME =======================================================================

    r_mod = [5094.02901670; 6127.87093630; 6380.24788850]

    # -- DCM -------------------------------------------------------------------------------

    D_teme_mod = r_eci_to_eci(MOD(), JD_UTC, TEME(), JD_UTC)

    r_teme = D_teme_mod * r_mod

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-7
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-7
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-7

    # -- Quaternion ------------------------------------------------------------------------

    q_teme_mod = r_eci_to_eci(Quaternion, MOD(), JD_UTC, TEME(), JD_UTC)

    r_teme = vect(q_teme_mod \ r_mod * q_teme_mod)

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-7
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-7
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-7

    # == TEME => MOD =======================================================================

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_mod_teme = r_eci_to_eci(TEME(), JD_UTC, MOD(), JD_UTC)

    r_mod = D_mod_teme * r_teme

    @test r_mod[1] ≈ +5094.02901670 atol = 1e-7
    @test r_mod[2] ≈ +6127.87093630 atol = 1e-7
    @test r_mod[3] ≈ +6380.24788850 atol = 1e-7

    # -- Quaternion ------------------------------------------------------------------------

    q_mod_teme = r_eci_to_eci(Quaternion, TEME(), JD_UTC, MOD(), JD_UTC)

    r_mod = vect(q_mod_teme \ r_teme * q_mod_teme)

    @test r_mod[1] ≈ +5094.02901670 atol = 1e-7
    @test r_mod[2] ≈ +6127.87093630 atol = 1e-7
    @test r_mod[3] ≈ +6380.24788850 atol = 1e-7
end

# == TOD <=> TEME ==========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#
# one gets:
#
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci TOD <=> TEME" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == TOD => TEME =======================================================================

    r_tod = [5094.51620300; 6127.36527840; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_teme_tod = r_eci_to_eci(TOD(), JD_UTC, TEME(), JD_UTC, eop_iau1980)

    r_teme = D_teme_tod * r_tod

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-4
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-4
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_teme_tod = r_eci_to_eci(Quaternion, TOD(), JD_UTC, TEME(), JD_UTC, eop_iau1980)

    r_teme = vect(q_teme_tod \ r_tod * q_teme_tod)

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-4
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-4
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-4

    # == TEME => TOD =======================================================================

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_tod_teme = r_eci_to_eci(TEME(), JD_UTC, TOD(), JD_UTC, eop_iau1980)

    r_tod = D_tod_teme * r_teme

    @test r_tod[1] ≈ +5094.51620300 atol = 1e-4
    @test r_tod[2] ≈ +6127.36527840 atol = 1e-4
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_tod_teme = r_eci_to_eci(Quaternion, TEME(), JD_UTC, TOD(), JD_UTC, eop_iau1980)

    r_tod = vect(q_tod_teme \ r_teme * q_tod_teme)

    @test r_tod[1] ≈ +5094.51620300 atol = 1e-4
    @test r_tod[2] ≈ +6127.36527840 atol = 1e-4
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-4

    # == No EOP corrections ================================================================

    # == TOD => TEME =======================================================================

    r_tod = [5094.51478040; 6127.36646120; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_teme_tod = r_eci_to_eci(TOD(), JD_UTC, TEME(), JD_UTC)

    r_teme = D_teme_tod * r_tod

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-4
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-4
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_teme_tod = r_eci_to_eci(Quaternion, TOD(), JD_UTC, TEME(), JD_UTC)

    r_teme = vect(q_teme_tod \ r_tod * q_teme_tod)

    @test r_teme[1] ≈ +5094.18016210 atol = 1e-4
    @test r_teme[2] ≈ +6127.64465950 atol = 1e-4
    @test r_teme[3] ≈ +6380.34453270 atol = 1e-4

    # == TEME => TOD =======================================================================

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_tod_teme = r_eci_to_eci(TEME(), JD_UTC, TOD(), JD_UTC)

    r_tod = D_tod_teme * r_teme

    @test r_tod[1] ≈ +5094.51478040 atol = 1e-7
    @test r_tod[2] ≈ +6127.36646120 atol = 1e-7
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-7

    # -- Quaternion ------------------------------------------------------------------------

    q_tod_teme = r_eci_to_eci(Quaternion, TEME(), JD_UTC, TOD(), JD_UTC)

    r_tod = vect(q_tod_teme \ r_teme * q_tod_teme)

    @test r_tod[1] ≈ +5094.51478040 atol = 1e-7
    @test r_tod[2] ≈ +6127.36646120 atol = 1e-7
    @test r_tod[3] ≈ +6380.34453270 atol = 1e-7
end

############################################################################################
#                                IAU-2006 / 2010 CIO-based                                 #
############################################################################################

# == GCRF <=> CIRS =========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
# NOTE: It seems that the results in the Example 3-14 is computed **without** the `dX` and
# `dY` corrections, whereas in the Table 3-6 they are computed **with** the corrections.
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_cirs = 5100.01840470  i + 6122.78636480 j + 6380.34453270 k [km]
#
# one gets the following (this is the result in Table 3-6):
#
#   r_gcrf = 5102.50895290  i + 6123.01139910 j + 6378.13693380 k [km]
#
############################################################################################

@testset "Function r_eci_to_eci GCRF <=> CIRS" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == GCRF => CIRS ======================================================================

    r_gcrf = [5102.50895290; 6123.01139910; 6378.13693380]

    # -- DCM -------------------------------------------------------------------------------

    D_cirs_gcrf = r_eci_to_eci(GCRF(), CIRS(), JD_UTC, eop_iau2000a)

    r_cirs = D_cirs_gcrf * r_gcrf

    @test r_cirs[1] ≈ +5100.01840470 atol = 1e-4
    @test r_cirs[2] ≈ +6122.78636480 atol = 1e-4
    @test r_cirs[3] ≈ +6380.34453270 atol = 1e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    D_cirs_gcrf = r_eci_to_eci(GCRF(), CIRS(), JD_UTC)

    r_cirs = D_cirs_gcrf * r_gcrf

    @test r_cirs[1] ≈ +5100.01840470 atol = 1e-4
    @test r_cirs[2] ≈ +6122.78636480 atol = 1e-4
    @test r_cirs[3] ≈ +6380.34453270 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_cirs_gcrf = r_eci_to_eci(Quaternion, GCRF(), CIRS(), JD_UTC, eop_iau2000a)

    r_cirs = vect(q_cirs_gcrf \ r_gcrf * q_cirs_gcrf)

    @test r_cirs[1] ≈ +5100.01840470 atol = 1e-4
    @test r_cirs[2] ≈ +6122.78636480 atol = 1e-4
    @test r_cirs[3] ≈ +6380.34453270 atol = 1e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    q_cirs_gcrf = r_eci_to_eci(Quaternion, GCRF(), CIRS(), JD_UTC)

    r_cirs = vect(q_cirs_gcrf \ r_gcrf * q_cirs_gcrf)

    @test r_cirs[1] ≈ +5100.01840470 atol = 1e-4
    @test r_cirs[2] ≈ +6122.78636480 atol = 1e-4
    @test r_cirs[3] ≈ +6380.34453270 atol = 1e-4

    # == CIRS => GCRF ======================================================================

    r_cirs = [+5100.01840470; +6122.78636480; +6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_gcrf_cirs = r_eci_to_eci(CIRS(), GCRF(), JD_UTC, eop_iau2000a)

    r_gcrf = D_gcrf_cirs * r_cirs

    @test r_gcrf[1] ≈ +5102.50895290 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01139910 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13693380 atol = 1e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    D_gcrf_cirs = r_eci_to_eci(CIRS(), GCRF(), JD_UTC)

    r_gcrf = D_gcrf_cirs * r_cirs

    @test r_gcrf[1] ≈ +5102.50895290 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01139910 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13693380 atol = 1e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_gcrf_cirs = r_eci_to_eci(Quaternion, CIRS(), GCRF(), JD_UTC, eop_iau2000a)

    r_gcrf = vect(q_gcrf_cirs \ r_cirs * q_gcrf_cirs)

    @test r_gcrf[1] ≈ +5102.50895290 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01139910 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13693380 atol = 1e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    q_gcrf_cirs = r_eci_to_eci(Quaternion, CIRS(), GCRF(), JD_UTC)

    r_gcrf = vect(q_gcrf_cirs \ r_cirs * q_gcrf_cirs)

    @test r_gcrf[1] ≈ +5102.50895290 atol = 1e-4
    @test r_gcrf[2] ≈ +6123.01139910 atol = 1e-4
    @test r_gcrf[3] ≈ +6378.13693380 atol = 1e-4
end

# == CIRS <=> CIRS =========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# A conversion between CIRS at the same epoch must return the identity rotation.
#
############################################################################################

@testset "Function r_eci_to_eci CIRS <=> CIRS" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == DCM ===============================================================================

    D_cirs_cirs = r_eci_to_eci(CIRS(), JD_UTC, CIRS(), JD_UTC)

    @test D_cirs_cirs[1, 1] ≈ 1 atol = 1e-14
    @test D_cirs_cirs[1, 2] ≈ 0 atol = 1e-14
    @test D_cirs_cirs[1, 3] ≈ 0 atol = 1e-14

    @test D_cirs_cirs[2, 1] ≈ 0 atol = 1e-14
    @test D_cirs_cirs[2, 2] ≈ 1 atol = 1e-14
    @test D_cirs_cirs[2, 3] ≈ 0 atol = 1e-14

    @test D_cirs_cirs[3, 1] ≈ 0 atol = 1e-14
    @test D_cirs_cirs[3, 2] ≈ 0 atol = 1e-14
    @test D_cirs_cirs[3, 3] ≈ 1 atol = 1e-14

    D_cirs_cirs = r_eci_to_eci(CIRS(), JD_UTC, CIRS(), JD_UTC, eop_iau2000a)

    @test D_cirs_cirs[1, 1] ≈ 1 atol = 1e-14
    @test D_cirs_cirs[1, 2] ≈ 0 atol = 1e-14
    @test D_cirs_cirs[1, 3] ≈ 0 atol = 1e-14

    @test D_cirs_cirs[2, 1] ≈ 0 atol = 1e-14
    @test D_cirs_cirs[2, 2] ≈ 1 atol = 1e-14
    @test D_cirs_cirs[2, 3] ≈ 0 atol = 1e-14

    @test D_cirs_cirs[3, 1] ≈ 0 atol = 1e-14
    @test D_cirs_cirs[3, 2] ≈ 0 atol = 1e-14
    @test D_cirs_cirs[3, 3] ≈ 1 atol = 1e-14

    # == Quaternion ========================================================================

    q_cirs_cirs = r_eci_to_eci(Quaternion, CIRS(), JD_UTC, CIRS(), JD_UTC)

    @test q_cirs_cirs[1] ≈ 1 atol = 1e-14
    @test q_cirs_cirs[2] ≈ 0 atol = 1e-14
    @test q_cirs_cirs[3] ≈ 0 atol = 1e-14
    @test q_cirs_cirs[4] ≈ 0 atol = 1e-14

    q_cirs_cirs = r_eci_to_eci(Quaternion, CIRS(), JD_UTC, CIRS(), JD_UTC, eop_iau2000a)

    @test q_cirs_cirs[1] ≈ 1 atol = 1e-14
    @test q_cirs_cirs[2] ≈ 0 atol = 1e-14
    @test q_cirs_cirs[3] ≈ 0 atol = 1e-14
    @test q_cirs_cirs[4] ≈ 0 atol = 1e-14
end

############################################################################################
#                              IAU-2006 / 2010 Equinox-based                               #
############################################################################################

# == GCRF <=> MJ2000 =======================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC     = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#
# one gets the following (this is the result in Table 3-6):
#
#   r_gcrf = 5102.50895780  i + 6123.01140380 j + 6378.13692530 k [km]
#
# Notice that the transformation we are testing here does not convert to the original J2000.
# However, this result is close enough for a test comparison.
#
############################################################################################

@testset "Function r_eci_to_eci GCRF <=> MJ2000" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == MJ2000 => GCRF ====================================================================

    r_mj2000 = [5102.50960000; 6123.01152000; 6378.13630000]

    # -- DCM -------------------------------------------------------------------------------

    D_GCRF_MJ2000 = r_eci_to_eci(MJ2000(), GCRF(), JD_UTC, eop_iau2000a)

    r_gcrf = D_GCRF_MJ2000 * r_mj2000

    @test r_gcrf[1] ≈ +5102.50895780 atol = 8e-4
    @test r_gcrf[2] ≈ +6123.01140380 atol = 8e-4
    @test r_gcrf[3] ≈ +6378.13692530 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    D_GCRF_MJ2000 = r_eci_to_eci(MJ2000(), GCRF(), JD_UTC)

    r_gcrf = D_GCRF_MJ2000 * r_mj2000

    @test r_gcrf[1] ≈ +5102.50895780 atol = 8e-4
    @test r_gcrf[2] ≈ +6123.01140380 atol = 8e-4
    @test r_gcrf[3] ≈ +6378.13692530 atol = 8e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_gcrf_mj2000 = r_eci_to_eci(Quaternion, MJ2000(), GCRF(), JD_UTC, eop_iau2000a)

    r_gcrf = vect(q_gcrf_mj2000 \ r_mj2000 * q_gcrf_mj2000)

    @test r_gcrf[1] ≈ +5102.50895780 atol = 8e-4
    @test r_gcrf[2] ≈ +6123.01140380 atol = 8e-4
    @test r_gcrf[3] ≈ +6378.13692530 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    q_gcrf_mj2000 = r_eci_to_eci(Quaternion, MJ2000(), GCRF(), JD_UTC)

    r_gcrf = vect(q_gcrf_mj2000 \ r_mj2000 * q_gcrf_mj2000)

    @test r_gcrf[1] ≈ +5102.50895780 atol = 8e-4
    @test r_gcrf[2] ≈ +6123.01140380 atol = 8e-4
    @test r_gcrf[3] ≈ +6378.13692530 atol = 8e-4

    # == GCRF => MJ2000 ====================================================================

    r_gcrf = [+5102.50895780; +6123.01140380; +6378.13692530]

    # -- DCM -------------------------------------------------------------------------------

    D_mj2000_gcrf = r_eci_to_eci(GCRF(), MJ2000(), JD_UTC, eop_iau2000a)

    r_mj2000 = D_mj2000_gcrf * r_gcrf

    @test r_mj2000[1] ≈ +5102.50960000 atol = 8e-4
    @test r_mj2000[2] ≈ +6123.01152000 atol = 8e-4
    @test r_mj2000[3] ≈ +6378.13630000 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    D_mj2000_gcrf = r_eci_to_eci(GCRF(), MJ2000(), JD_UTC)

    r_mj2000 = D_mj2000_gcrf * r_gcrf

    @test r_mj2000[1] ≈ +5102.50960000 atol = 8e-4
    @test r_mj2000[2] ≈ +6123.01152000 atol = 8e-4
    @test r_mj2000[3] ≈ +6378.13630000 atol = 8e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_mj2000_gcrf = r_eci_to_eci(Quaternion, GCRF(), MJ2000(), JD_UTC, eop_iau2000a)

    r_mj2000 = vect(q_mj2000_gcrf \ r_gcrf * q_mj2000_gcrf)

    @test r_mj2000[1] ≈ +5102.50960000 atol = 8e-4
    @test r_mj2000[2] ≈ +6123.01152000 atol = 8e-4
    @test r_mj2000[3] ≈ +6378.13630000 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    q_mj2000_gcrf = r_eci_to_eci(Quaternion, GCRF(), MJ2000(), JD_UTC)

    r_mj2000 = vect(q_mj2000_gcrf \ r_gcrf * q_mj2000_gcrf)

    @test r_mj2000[1] ≈ +5102.50960000 atol = 8e-4
    @test r_mj2000[2] ≈ +6123.01152000 atol = 8e-4
    @test r_mj2000[3] ≈ +6378.13630000 atol = 8e-4
end

# == GCRF <=> MOD ==========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895780  i + 6123.01140380 j + 6378.13692530 k [km]
#
# one gets the following (this is the result in Table 3-6):
#
#   r_mod  = +5094.02896110   i + 6127.87113500   j + 6380.24774200   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci GCRF <=> MOD" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == MOD => GCRF =======================================================================

    r_mod = [+5094.02896110; +6127.87113500; +6380.24774200]

    # -- DCM -------------------------------------------------------------------------------

    D_gcrf_mod = r_eci_to_eci(MOD06(), GCRF(), JD_UTC, eop_iau2000a)

    r_gcrf = D_gcrf_mod * r_mod

    @test r_gcrf[1] ≈ +5102.50895780 atol = 1e-7
    @test r_gcrf[2] ≈ +6123.01140380 atol = 1e-7
    @test r_gcrf[3] ≈ +6378.13692530 atol = 1e-7

    # -- Quaternion ------------------------------------------------------------------------

    q_gcrf_mod = r_eci_to_eci(Quaternion, MOD06(), GCRF(), JD_UTC, eop_iau2000a)

    r_gcrf = vect(q_gcrf_mod \ r_mod * q_gcrf_mod)

    @test r_gcrf[1] ≈ +5102.50895780 atol = 1e-7
    @test r_gcrf[2] ≈ +6123.01140380 atol = 1e-7
    @test r_gcrf[3] ≈ +6378.13692530 atol = 1e-7

    # == GCRF => MOD =======================================================================

    r_gcrf = [+5102.50895780; +6123.01140380; +6378.13692530]

    # -- DCM -------------------------------------------------------------------------------

    D_mod_gcrf = r_eci_to_eci(GCRF(), MOD06(), JD_UTC, eop_iau2000a)

    r_mod = D_mod_gcrf * r_gcrf

    @test r_mod[1] ≈ +5094.02896110 atol = 1e-7
    @test r_mod[2] ≈ +6127.87113500 atol = 1e-7
    @test r_mod[3] ≈ +6380.24774200 atol = 1e-7

    # -- Quaternion ------------------------------------------------------------------------

    q_mod_gcrf = r_eci_to_eci(Quaternion, GCRF(), MOD06(), JD_UTC, eop_iau2000a)

    r_mod = vect(q_mod_gcrf \ r_gcrf * q_mod_gcrf)

    @test r_mod[1] ≈ +5094.02896110 atol = 1e-7
    @test r_mod[2] ≈ +6127.87113500 atol = 1e-7
    @test r_mod[3] ≈ +6380.24774200 atol = 1e-7
end

# == GCRF <=> ERS ==========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895780  i + 6123.01140380 j + 6378.13692530 k [km]
#
# one gets the following (this is the result in Table 3-6):
#
#   r_ers  = +5094.51462800   i + 6127.36658790   j + 6380.34453270   k [km]
#
############################################################################################

@testset "Function r_eci_to_eci GCRF <=> ERS" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == ERS => GCRF =======================================================================

    r_ers = [+5094.51462800; +6127.36658790; +6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_gcrf_ers = r_eci_to_eci(ERS(), GCRF(), JD_UTC, eop_iau2000a)

    r_gcrf = D_gcrf_ers * r_ers

    @test r_gcrf[1] ≈ +5102.50895780 atol = 5e-5
    @test r_gcrf[2] ≈ +6123.01140380 atol = 5e-5
    @test r_gcrf[3] ≈ +6378.13692530 atol = 5e-5

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    D_gcrf_ers = r_eci_to_eci(ERS(), GCRF(), JD_UTC)

    r_gcrf = D_gcrf_ers * r_ers

    @test r_gcrf[1] ≈ +5102.50895780 atol = 5e-5
    @test r_gcrf[2] ≈ +6123.01140380 atol = 5e-5
    @test r_gcrf[3] ≈ +6378.13692530 atol = 5e-5

    # -- Quaternion ------------------------------------------------------------------------

    q_gcrf_ers = r_eci_to_eci(Quaternion, ERS(), GCRF(), JD_UTC, eop_iau2000a)

    r_gcrf = vect(q_gcrf_ers \ r_ers * q_gcrf_ers)

    @test r_gcrf[1] ≈ +5102.50895780 atol = 5e-5
    @test r_gcrf[2] ≈ +6123.01140380 atol = 5e-5
    @test r_gcrf[3] ≈ +6378.13692530 atol = 5e-5

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    q_gcrf_ers = r_eci_to_eci(Quaternion, ERS(), GCRF(), JD_UTC)

    r_gcrf = vect(q_gcrf_ers \ r_ers * q_gcrf_ers)

    @test r_gcrf[1] ≈ +5102.50895780 atol = 5e-5
    @test r_gcrf[2] ≈ +6123.01140380 atol = 5e-5
    @test r_gcrf[3] ≈ +6378.13692530 atol = 5e-5

    # == GCRF => ERS =======================================================================

    r_gcrf = [+5102.50895780; +6123.01140380; +6378.13692530]

    # -- DCM -------------------------------------------------------------------------------

    D_ers_gcrf = r_eci_to_eci(GCRF(), ERS(), JD_UTC, eop_iau2000a)

    r_ers = D_ers_gcrf * r_gcrf

    @test r_ers[1] ≈ +5094.51462800 atol = 5e-5
    @test r_ers[2] ≈ +6127.36658790 atol = 5e-5
    @test r_ers[3] ≈ +6380.34453270 atol = 5e-5

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    D_ers_gcrf = r_eci_to_eci(GCRF(), ERS(), JD_UTC)

    r_ers = D_ers_gcrf * r_gcrf

    @test r_ers[1] ≈ +5094.51462800 atol = 5e-5
    @test r_ers[2] ≈ +6127.36658790 atol = 5e-5
    @test r_ers[3] ≈ +6380.34453270 atol = 5e-5

    # -- Quaternion ------------------------------------------------------------------------

    q_ers_gcrf = r_eci_to_eci(Quaternion, GCRF(), ERS(), JD_UTC, eop_iau2000a)

    r_ers = vect(q_ers_gcrf \ r_gcrf * q_ers_gcrf)

    @test r_ers[1] ≈ +5094.51462800 atol = 5e-5
    @test r_ers[2] ≈ +6127.36658790 atol = 5e-5
    @test r_ers[3] ≈ +6380.34453270 atol = 5e-5

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    q_ers_gcrf = r_eci_to_eci(Quaternion, GCRF(), ERS(), JD_UTC)

    r_ers = vect(q_ers_gcrf \ r_gcrf * q_ers_gcrf)

    @test r_ers[1] ≈ +5094.51462800 atol = 5e-5
    @test r_ers[2] ≈ +6127.36658790 atol = 5e-5
    @test r_ers[3] ≈ +6380.34453270 atol = 5e-5
end

# == MJ2000 <=> MOD ========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC     = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#
# one gets the following (this is the result in Table 3-6):
#
#   r_mod  = +5094.02896110   i + 6127.87113500   j + 6380.24774200   k [km]
#
# Notice that the transformation we are testing here does not convert to the original J2000.
# However, this result is close enough for a test comparison.
#
############################################################################################

@testset "Function r_eci_to_eci MJ2000 <=> MOD" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == MOD => MJ2000 =====================================================================

    r_mod = [+5094.02896110; +6127.87113500; +6380.24774200]

    # -- DCM -------------------------------------------------------------------------------

    D_mj2000_mod = r_eci_to_eci(MOD06(), MJ2000(), JD_UTC, eop_iau2000a)

    r_mj2000 = D_mj2000_mod * r_mod

    @test r_mj2000[1] ≈ +5102.50960000 atol = 8e-4
    @test r_mj2000[2] ≈ +6123.01152000 atol = 8e-4
    @test r_mj2000[3] ≈ +6378.13630000 atol = 8e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_mj2000_mod = r_eci_to_eci(Quaternion, MOD06(), MJ2000(), JD_UTC, eop_iau2000a)

    r_mj2000 = vect(q_mj2000_mod \ r_mod * q_mj2000_mod)

    @test r_mj2000[1] ≈ +5102.50960000 atol = 8e-4
    @test r_mj2000[2] ≈ +6123.01152000 atol = 8e-4
    @test r_mj2000[3] ≈ +6378.13630000 atol = 8e-4

    # == MJ2000 => MOD =====================================================================

    r_mj2000 = [5102.50960000; 6123.01152000; 6378.13630000]

    # -- DCM -------------------------------------------------------------------------------

    D_MOD_MJ2000 = r_eci_to_eci(MJ2000(), MOD06(), JD_UTC, eop_iau2000a)

    r_mod = D_MOD_MJ2000 * r_mj2000

    @test r_mod[1] ≈ +5094.02896110 atol = 8e-4
    @test r_mod[2] ≈ +6127.87113500 atol = 8e-4
    @test r_mod[3] ≈ +6380.24774200 atol = 8e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_mod_mj2000 = r_eci_to_eci(Quaternion, MJ2000(), MOD06(), JD_UTC, eop_iau2000a)

    r_mod = vect(q_mod_mj2000 \ r_mj2000 * q_mod_mj2000)

    @test r_mod[1] ≈ +5094.02896110 atol = 8e-4
    @test r_mod[2] ≈ +6127.87113500 atol = 8e-4
    @test r_mod[3] ≈ +6380.24774200 atol = 8e-4
end

# == MJ2000 <=> ERS ========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC     = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#
# one gets the following (this is the result in Table 3-6):
#
#   r_ers  = +5094.51462800   i + 6127.36658790   j + 6380.34453270   k [km]
#
# Notice that the transformation we are testing here does not convert to the original J2000.
# However, this result is close enough for a test comparison.
#
############################################################################################

@testset "Function r_eci_to_eci MJ2000 <=> ERS" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == ERS => MJ2000 =====================================================================

    r_ers = [+5094.51462800; +6127.36658790; +6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_mj2000_ers = r_eci_to_eci(ERS(), MJ2000(), JD_UTC, eop_iau2000a)

    r_mj2000 = D_mj2000_ers * r_ers

    @test r_mj2000[1] ≈ +5102.50960000 atol = 8e-4
    @test r_mj2000[2] ≈ +6123.01152000 atol = 8e-4
    @test r_mj2000[3] ≈ +6378.13630000 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    D_mj2000_ers = r_eci_to_eci(ERS(), MJ2000(), JD_UTC)

    r_mj2000 = D_mj2000_ers * r_ers

    @test r_mj2000[1] ≈ +5102.50960000 atol = 8e-4
    @test r_mj2000[2] ≈ +6123.01152000 atol = 8e-4
    @test r_mj2000[3] ≈ +6378.13630000 atol = 8e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_mj2000_ers = r_eci_to_eci(Quaternion, ERS(), MJ2000(), JD_UTC, eop_iau2000a)

    r_mj2000 = vect(q_mj2000_ers \ r_ers * q_mj2000_ers)

    @test r_mj2000[1] ≈ +5102.50960000 atol = 8e-4
    @test r_mj2000[2] ≈ +6123.01152000 atol = 8e-4
    @test r_mj2000[3] ≈ +6378.13630000 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    q_mj2000_ers = r_eci_to_eci(Quaternion, ERS(), MJ2000(), JD_UTC)

    r_mj2000 = vect(q_mj2000_ers \ r_ers * q_mj2000_ers)

    @test r_mj2000[1] ≈ +5102.50960000 atol = 8e-4
    @test r_mj2000[2] ≈ +6123.01152000 atol = 8e-4
    @test r_mj2000[3] ≈ +6378.13630000 atol = 8e-4

    # == MJ2000 => ERS =====================================================================

    r_mj2000 = [5102.50960000; 6123.01152000; 6378.13630000]

    # -- DCM -------------------------------------------------------------------------------

    D_ERS_MJ2000 = r_eci_to_eci(MJ2000(), ERS(), JD_UTC, eop_iau2000a)

    r_ers = D_ERS_MJ2000 * r_mj2000

    @test r_ers[1] ≈ +5094.51462800 atol = 8e-4
    @test r_ers[2] ≈ +6127.36658790 atol = 8e-4
    @test r_ers[3] ≈ +6380.34453270 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    D_ERS_MJ2000 = r_eci_to_eci(MJ2000(), ERS(), JD_UTC)

    r_ers = D_ERS_MJ2000 * r_mj2000

    @test r_ers[1] ≈ +5094.51462800 atol = 8e-4
    @test r_ers[2] ≈ +6127.36658790 atol = 8e-4
    @test r_ers[3] ≈ +6380.34453270 atol = 8e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_ers_mj2000 = r_eci_to_eci(Quaternion, MJ2000(), ERS(), JD_UTC, eop_iau2000a)

    r_ers = vect(q_ers_mj2000 \ r_mj2000 * q_ers_mj2000)

    @test r_ers[1] ≈ +5094.51462800 atol = 8e-4
    @test r_ers[2] ≈ +6127.36658790 atol = 8e-4
    @test r_ers[3] ≈ +6380.34453270 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    q_ers_mj2000 = r_eci_to_eci(Quaternion, MJ2000(), ERS(), JD_UTC)

    r_ers = vect(q_ers_mj2000 \ r_mj2000 * q_ers_mj2000)

    @test r_ers[1] ≈ +5094.51462800 atol = 8e-4
    @test r_ers[2] ≈ +6127.36658790 atol = 8e-4
    @test r_ers[3] ≈ +6380.34453270 atol = 8e-4
end

# == MOD <=> ERS ===========================================================================

############################################################################################
#                                       Test Results                                       #
############################################################################################
#
# == Scenario 01 ===========================================================================
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC     = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#
# one gets the following (this is the result in Table 3-6):
#
#   r_ers  = +5094.51462800   i + 6127.36658790   j + 6380.34453270   k [km]
#
# Notice that the transformation we are testing here does not convert to the original J2000.
# However, this result is close enough for a test comparison.
#
############################################################################################

@testset "Function r_eci_to_eci MOD <=> ERS" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    # == ERS => MOD ========================================================================

    r_ers = [+5094.51462800; +6127.36658790; +6380.34453270]

    # -- DCM -------------------------------------------------------------------------------

    D_mod_ers = r_eci_to_eci(ERS(), JD_UTC, MOD06(), JD_UTC, eop_iau2000a)

    r_mod = D_mod_ers * r_ers

    @test r_mod[1] ≈ +5094.02896110 atol = 8e-4
    @test r_mod[2] ≈ +6127.87113500 atol = 8e-4
    @test r_mod[3] ≈ +6380.24774200 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    D_mod_ers = r_eci_to_eci(ERS(), JD_UTC, MOD06(), JD_UTC)

    r_mod = D_mod_ers * r_ers

    @test r_mod[1] ≈ +5094.02896110 atol = 8e-4
    @test r_mod[2] ≈ +6127.87113500 atol = 8e-4
    @test r_mod[3] ≈ +6380.24774200 atol = 8e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_mod_ers = r_eci_to_eci(Quaternion, ERS(), JD_UTC, MOD06(), JD_UTC, eop_iau2000a)

    r_mod = vect(q_mod_ers \ r_ers * q_mod_ers)

    @test r_mod[1] ≈ +5094.02896110 atol = 8e-4
    @test r_mod[2] ≈ +6127.87113500 atol = 8e-4
    @test r_mod[3] ≈ +6380.24774200 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    q_mod_ers = r_eci_to_eci(Quaternion, ERS(), JD_UTC, MOD06(), JD_UTC)

    r_mod = vect(q_mod_ers \ r_ers * q_mod_ers)

    @test r_mod[1] ≈ +5094.02896110 atol = 8e-4
    @test r_mod[2] ≈ +6127.87113500 atol = 8e-4
    @test r_mod[3] ≈ +6380.24774200 atol = 8e-4

    # == MOD => ERS ========================================================================

    r_mod = [+5094.02896110; +6127.87113500; +6380.24774200]

    # -- DCM -------------------------------------------------------------------------------

    D_ers_mod = r_eci_to_eci(MOD06(), JD_UTC, ERS(), JD_UTC, eop_iau2000a)

    r_ers = D_ers_mod * r_mod

    @test r_ers[1] ≈ +5094.51462800 atol = 8e-4
    @test r_ers[2] ≈ +6127.36658790 atol = 8e-4
    @test r_ers[3] ≈ +6380.34453270 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    D_ers_mod = r_eci_to_eci(MOD06(), JD_UTC, ERS(), JD_UTC)

    r_ers = D_ers_mod * r_mod

    @test r_ers[1] ≈ +5094.51462800 atol = 8e-4
    @test r_ers[2] ≈ +6127.36658790 atol = 8e-4
    @test r_ers[3] ≈ +6380.34453270 atol = 8e-4

    # -- Quaternion ------------------------------------------------------------------------

    q_ers_mod = r_eci_to_eci(Quaternion, MOD06(), JD_UTC, ERS(), JD_UTC, eop_iau2000a)

    r_ers = vect(q_ers_mod \ r_mod * q_ers_mod)

    @test r_ers[1] ≈ +5094.51462800 atol = 8e-4
    @test r_ers[2] ≈ +6127.36658790 atol = 8e-4
    @test r_ers[3] ≈ +6380.34453270 atol = 8e-4

    # NOTE: We do not have the values without the EOP corrections. However, the difference
    # is very small for this case and we will test using the same values.
    q_ers_mod = r_eci_to_eci(Quaternion, MOD06(), JD_UTC, ERS(), JD_UTC)

    r_ers = vect(q_ers_mod \ r_mod * q_ers_mod)

    @test r_ers[1] ≈ +5094.51462800 atol = 8e-4
    @test r_ers[2] ≈ +6127.36658790 atol = 8e-4
    @test r_ers[3] ≈ +6380.34453270 atol = 8e-4
end
